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Abstract 

We describe the software package FELIX that solves the equations of the time-dependent generator coordinate method (TDGCM) in 
A'-dimensions (N > 1) under the Gaussian overlap approximation. The numerical resolution is based on the Galerkin finite element 
discretization of the collective space and the Crank-Nicolson scheme for time integration. The TDGCM solver is implemented 
entirely in C++. Several additional tools written in C++, Python or bash scripting language are also included for convenience. In 
this paper, the solver is tested with a series of benchmarks calculations. We also demonstrate the ability of our code to handle a 
realistic calculation of fission dynamics. 
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PROGRAM SUMMARY/NEW VERSION PROGRAM 
SUMMARY 


Program Title: FELIX-1.0 
Journal Reference: 

Catalogue identifier: 

Licensing provisions: 

Programming language: C++ 

Computer: Intel Xeon, Intel Core 
Operating system: LINUX 

RAM: Memory usage depends on the number of nodes in the 
calculation mesh as well as on the degree of the interpolation 
polynomials. For a ID calculation with linear polynomials on a mesh 
with 600 nodes, memory usage is approximately 3.3 MB; in a realistic 
simulation of fission on a 2D mesh with quadratic polynomials and 
1.3 10 s nodes, it reaches 1.5 GiB. 

Number of processors used: 

The code is multi-threaded based on the OpenMP API specification 
for parallel programming. Any number of threads may be specified by 
the user. 

Keywords: FELIX; Finite element method; Generator coordinate 
method; Gaussian overlap approximation; Nuclear fission; 
Classification: 17.23 Fission and Fusion Processes 
External routines/libraries: The solver itself requires the BLAS and 
LAPACK libraries, and a Fortran compiler with OpenMP support. 
Building the documentation requires DoxyGen-1.8.6 or higher. Build¬ 
ing the full set of tools also requires GSL, PETSc, SLEPc and Boost. 
In particular, environment variables PETSC DHL PF.TSC ARCH, 
SLEPC_DIR and SLEPC_ARCH must be set. 

Nature of problem: 

Nuclear fission is a relatively slow process compared to the typical 
timescale of the intrinsic motion of the nucleons. In the adiabatic 
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approximation, it can be described as a large amplitude collective 
motion driven by only a few collective degrees of freedom. In the 
time-dependent generator coordinate method (TDGCM), the nuclear 
wave-function is thus described as a time-dependent, linear superpo¬ 
sition of basis functions in this collective space. Further assuming a 
Gaussian overlap approximation (GOA) for the basis functions, the 
time-dependent Schrodinger equation can be reduced into a local, 
time-dependent, Schrodinger-like equation in collective space. This 
is the TDGCM+GOA equation. Scission configurations are defined 
as a hyper-surface in the A-dimensional collective space. Fission 
fragment distributions are then computed by integrating over time the 
flux of the collective wave-packet across the scission hyper-surface. 
This microscopic approach to fission fragment distributions is fully 
quantum-mechanical. 

Solution method: 

FELIX solves the TDGCM+GOA equation by using the Galerkin 
finite element method to discretize the A-dimensional collective 
space, and the Crank-Nicolson scheme to solve for the time evolution. 
At each time step, this procedure requires solving a linear system 
of equation involving sparse, complex, symmetric matrices. FELIX 
employs an iterative QMR algorithm to perform matrix inversion. 
Restrictions: 

Although the program can operate in an arbitrary number of dimen¬ 
sions N, it has only been tested in practice on 1, 2 and 3 dimensional 
meshes. 

Unusual features: 

Additional comments: 

The code has checkpointing capabilities: the collective wave-function, 
norm a and energy kernels are stored on disk every n iterations, 
ensuring that the program can resume where it stops. 

Running time: 

Running time grows linearly with the number of time-steps requested 
by the user. It is also highly dependent on the number of nodes in the 
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space mesh. Two periods of a ID harmonic oscillator (600 nodes, 800 
time steps) are typically computed in a few seconds on one thread of a 
Intel(R) Core(TM) i5 CPU. A 2-dimensional realistic case of fission 
(10 s nodes, 10 5 time steps) requires roughly 10 hours on 10 threads of 
an Intel Xeon EP X5660 processor. 


1. Introduction 

Induced nuclear fission plays an essential role in important 
societal applications ranging from stockpile science to criti¬ 
cal assemblies for new generation nuclear reactors [1], It is 
also one of the leading mechanisms determining the stability 
of super-heavy elements and the end point of nucleosynthesis 
[2] . Many of these applications require the detailed knowledge 
of fission product yields (FPY), which may include the charge, 
mass, kinetic energy and excitation energy distribution of the 
fission fragments. In spite of recent technological advances, 
FPY measurements are not always possible, especially in very 
short-lived neutron-rich or heavy nuclei. Predictions based on 
theoretical models of fission are thus unavoidable. 

While there exists a number of powerful phenomenological 
or semi-microscopic models on the market, the consensus is 
that a truly predictive theory of fission should ultimately be 
based only on our knowledge of nuclear forces and quantum 
many-body methods. In this context, most of the effort in 
the last decades has been focused on describing the dynam¬ 
ics of induced fission in the time-dependent generator coor¬ 
dinate method (TDGCM) associated with the Gaussian over¬ 
lap approximation (GOA) [3, 4, 5]. Under these assumptions, 
the original time-dependent many-body Schrodinger equation 
is reduced to a local Schrodinger-like equation that depends on 
only a few relevant collective variables. This approach was 
able to predict the characteristic times of low-energy induced 
fission [6, 7]. More recently, it was also successfully used to 
provide the first estimate of the mass and kinetic energy distri¬ 
butions of fission fragments for neutron-induced fission in the 
actinides region [8, 9, 10]. 

Until now, the aforementioned calculations have been based 
on the discretization of the TDGCM equations using finite dif¬ 
ferences in a regularly meshed hyper-cube. Given the computa¬ 
tional resources required by this simple scheme, fission dynam¬ 
ics has only been studied in 2-dimensional collective spaces. 
Yet, it is well-known from both semi-phenomenological and 
fully microscopic approaches that at least four or five collective 
variables play a role in the dynamics of fission [11, 12, 13, 14]. 
Although possible in theory, extending the current scheme to 
N > 2 collective spaces would be prohibitive computation¬ 
ally. Similarly, increasing the fidelity of the calculation from 
no points to n points for all N dimensions scales approximately 
like ( n/no) N . A lot of this increase in computational cost would 
be wasted in regions of the collective space far away from the 
scission configurations, where the resolution does not need to 
be very high. For these reasons, it is highly desirable to move 
to a more flexible, more scalable discretization scheme. 

In this work, we thus introduce the code FELIX, which im¬ 
plements the Galerkin finite element method for the discretiza¬ 


tion of the TDGCM equation. This well-known method allows 
the use of irregularly-spaced meshes in the collective space, 
which results in turn in calculations scaling much more effi¬ 
ciently with the number of collective variables. In addition, 
standard p-refinement techniques, i.e., the use of higher-degree 
polynomial bases in each finite element, give a better control 
on the numerical precision of the calculations. Finally, there 
are virtually no restrictions in the number of collective variables 
used. The code FELIX has been tested up to N = 3. 

After a brief introduction to the TDGCM+GOA equations in 
section 2, we present in section 3 the numerical implementation 
in the code FELIX. The implementation is validated through a 
series of simple benchmarks discussed in section 4. The conver¬ 
gence of a realistic calculation of the fission of 240 Pu is demon¬ 
strated in section 5. Finally, sections 6 and 7 give practical 
information on how to install and use the code FELIX. 


2. Fission dynamics in the TDGCM+GOA approach 

In this section, we briefly recall how to obtain a collective, 
Schrodinger-like equation to describe low-energy nuclear dy¬ 
namics. We also explain how fission fragment distributions can 
be extracted from the integration of the collective flux across 
the hyper-surface defining the scission configurations. 


2.1. The TDGCM+GOA equation 

We recall that the time evolution of a many-body quantum 
system is given by the time-dependent Schrodinger equation, 
which is obtained from the variation <kS['P] = 0 of the quantum 
mechanical action given by [15] 


sm = 



mo i [H-nif] mo> 
momo> 


(i) 


where | V P(0) is the full many-body wave-function for the sys¬ 
tem. In most nuclear physics applications, the nuclear Hamil¬ 
tonian H contains an effective two-body potential such as, e.g., 
the Skyrme or Gogny interaction. In the time-dependent gen¬ 
erator coordinate method (TDGCM), the nuclear many-body 
wave function |'T(f)) takes the form[3, 4, 5] 

m?)> = fdq f{q,t)mq)). (2) 

The functions m^r)) are known many-body states parametrized 
by a vector of collective variables q. In the context of fission, 
the m^r)) are chosen as the solutions to the static Hartree-Fock- 
Bogoliubov (HFB) equations under a set of constraints q. These 
constraints, which are the collective variables driving the fis¬ 
sion process, can be expectation values of multipole moments, 
quantities related to pairing such as particle number fluctuations 
AN 2 , etc. Recall that the HFB solutions at point q are entirely 
characterized by the one-body density matrix p and two-body 
pairing tensor k. 

Inserting the ansatz of Eq.(2) in the variational principle (1) 
yields the (time-dependent) Hill-Wheeler equation. In contrast 
to the static case, there has been no attempt so far to solve 
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the time-dependent Hill-Wheeler equation numerically, as the 
computational resources needed are beyond current capabili¬ 
ties. Instead, a widespread approach consists in assuming that 
the norm kernels (TT^OITT#')) can be approximated by a Gaus¬ 
sian form factor [16]. Inserting this Gaussian overlap approxi¬ 
mation (GOA) into the Hill-Wheeler equation (using a second 
order expansion in q - q ') leads to a local, time-dependent, 
Schrodinger-like equation in the space Q of collective coordi¬ 
nates q. 


d 


tf d d 

— V + v(q) 

2 dq k dqi 


g(q,t), (3) 


where 

• The function g(q, t ) is complex. It is related to the weight 
function f(q, t ) appearing in Eq.(2) and contains all the 
information about the dynamics of the system. Moreover, 
the quantity | g(q, t)\ 2 can be interpreted as the probability 
density for the system to be in the state |T , (</)) at time f; 
see also section 2.2 below. 

• The real scalar field V(q) and the real symmetric tensor 
field Bkiiq) are fully determined by the knowledge of the 
effective Hamiltonian H and the generator states \'V(q)}. 
They reflect the static nuclear properties of the system un¬ 
der study. 


- the external region. The hyper-surface separating the two re¬ 
gions corresponds to the set of scission configurations. The rig¬ 
orous definition and accurate determination of these scission 
configurations are themselves challenging problems, which go 
beyond the scope of this paper, see Refs. [17, 18, 19, 12, 20, 14] 
for additional discussions. For practical calculations of fission 
fragment distributions with FELIX we will simply assume the 
existence of such a scission hyper-surface. 

In general, the local, one-body density matrix p{r) in each of 
the scission points in the collective space Q is characterized by 
two high-density regions separated by a thin neck. Assuming 
the neck is located along the z-axis of the intrinsic reference 
frame, the charge and mass of each fragment can be obtained 
by simple integration of p(r) over the domains z e] - oo, z,v] 
and z e [zjv, +°°[; see, e.g., [19, 12, 14], According to this 
procedure, one can associate with each point q of the scission 
hyper-surface a pair of fragment masses. It follows that the 
flux of the probability current (6) through the scission hyper¬ 
surface gives a very good estimate of the relative probability of 
observing a given pair of fragments at time t. We thus define 
the integrated flux F(£, 1) through an oriented surface element £ 
as ^ 

F(£,t) = f dT f J(q, t) ■ dS. (7) 

Jt =0 Jqzt; 

Following [7, 10], the fission fragment mass yield for mass A is 
defined formally as 


Throughout this paper, equation (3) will be referred to as the 
TDGCM+GOA equation. 


Y(A ) oc V lim F {£,£), 

I J t —> + oo 




( 8 ) 


2.2. Collective flux and fission fragment distributions 

The TDGCM+GOA equation implies a continuity equation 
for the probability density | g(q, t )| 2 , 

11 g(q,t)\ 2 = -V-J(q,t). (4) 

at 

The real vector field J(q, t ) is thus a current of probability. It 
can be expressed formally as a function of the collective wave- 
function, 

J(q, t ) = £ B(q ) [g\q, t)Wg(q, t) - g(q , t)Vg\q, f)]. (5) 

2i 

Specifically, the coordinates of the current of probability read 


Jk(q,t) = 7T. y Bkfq) 
2l tt 


g*(q-t)^-(q,t) - g(q,t)^-(q,t) 
oqi dqi 


( 6 ) 

As our system evolves in time, its density probability will flow 
starting from the area of the collective space Q where the initial 
wave-function was localized. This evolution is driven by the 
Hamiltonian H through the inertia tensor Bifq) and the poten¬ 
tial energy surface V(q). 

In the case of fission, the potential energy surface is com¬ 
puted up to the points q where the nuclear geometry corre¬ 
sponds to two well-separated fragments. One can thus partition 
the space Q into a region where the nucleus is whole - the in¬ 
ternal region, and another where it has split in two fragments 


where F\ is the set of all oriented hyper-surfaces f belonging 
to the scission hyper-surface such that one of the fragments has 
mass A. In practice, our calculation of the fragments mass num¬ 
ber produces non integer values. Moreover, one elementary 
surface £ may contain several fragmentations. In this work, 
we equally distribute the flux component F(£, t) between the 
masses calculated at the vertices of the edge £: 

Y( A ) = C V -J- V lim F(£,t), (9) 

Z —l N t—l r—>+oo 
f ve.7l(f) 

The sum on £ runs on the whole scission hyper-surface. The set 
£A(£) contains the vertices of £ at which one of the fragments 
has a mass in the interval [A - 1 /2; A + 1 /2\. The normalization 
constant C is chosen as usual such that 

A total 

y Y(A ) = 200. (10) 

A=0 

In practice, the flux is only integrated from t = 0 to t = f max . 
Equations (6)-(10) show how to extract fission fragment yields 
from the knowledge of the collective wave function g(q, t) so¬ 
lution to the TDGCM+GOA equations. 

3. Numerical methods 

In this section, we detail the numerical methods implemented 
in the code FELIX to solve the TDGCM+GOA equation (3) and 
calculate the flux defined in Eq. (7). 
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3.1. Restriction to a finite domain of space 

FELIX solves Eq. (3) in a finite domain Q of the collective 
space Q. To ensure the uniqueness of the solution, Dirichlet 
conditions are imposed at the boundary dD. of the domain, 


With this definition, Eq. (12) can be recast into 

Vf> e £ 2 (Q C), Vf e [0, f max ] : <*K0> = 0, (15) 

with the residual r(q, t) defined as 


Vqedd: g(q,t) = 0. (11) 

Imposing this condition is justified as long as the actual solu¬ 
tion g(q.t) is well confined inside the domain Q during the 
whole time evolution of the system. In practice, this may re¬ 
quire choosing an excessively large domain £1 In the case of 
fission for example, only the internal region discussed in sec¬ 
tion 2.2 and its interface with the external zone present a phys¬ 
ical interest. However, we cannot limit Q. to this area because 
the probability to observe the fissioning system outside of this 
configuration subset is not negligible. 

To circumvent this issue, FELIX defines an absorption band 
along the boundary <9£2. This band artificially simulates the 
leakage of the wave packet g(q, t) outside of the calculation do¬ 
main. Formally, absorption is taken into account by introducing 
a new imaginary term in the evolution equation, 


V q £ Q, t £ [0, f max ] : 




12 Q Q 

-y ^ dq~ k Bkl ^d^, + V(q ^ ~ iM(?) 


dq, 


g(q, t). 
( 12 ) 


The real scalar field A(q ) is non zero only in the absorption 
band. In this region, A(q) is taken as a simple polynomial in¬ 
creasing smoothly from 0 on the inner border of the band and 
reaching its maximum at the boundary of the domain. 


A(q) = 4 r 



(13) 


r(q.t) 


h 2 v—i d d 

2 t-f dq k dqi 


-ihA(q) - ih — 
ot 


g(q,t) (16) 


Following the standard approach of the finite elements 
method, the domain Q is first partitioned into a mesh. In our 
case, each cell of the mesh is a IV-dimensional simplex (tri¬ 
angle if N = 2, tetrahedron if A' = 3, etc.). We note 5 the 
set of all simplices in the domain. Inside every simplex of the 
mesh, we assume a polynomial form for the numerical solution 
of Eq. (12). At any time t and in any simplex s £ S. we thus 
define the local interpolating polynomial P s , 


Vs £ S,Vq £ s : g(q,t) = P s ,,(q). (17) 


For each simplex ,v e S, we select the degree d s of the interpo¬ 
lating polynomial. The space *P S of all interpolating polynomi¬ 
als in the simplex s € 5 is a vector space. Its dimension D s is 
given by the binomial coefficient, 


Os = 


N + d s \ 
d s )■ 


(18) 


In order to discretize Eq. (12), we now build a convenient basis 
of the space V s . First, we define for each simplex s £ S a finite 
set of specific points y £ s called nodes. Next, we introduce a 
set of real polynomials </> s j associated with the simplex s. For 
all nodes i of the simplex s, the polynomial cf> x j is defined by 
the requirement 


The quantity x(q) is the minimal Euclidean distance between 
the point q and the boundary dCl. The parameters r and w cor¬ 
respond to the average absorption rate and width of the absorp¬ 
tion band respectively. These two parameters can be tuned by 
the user as a function of the problem characteristics to ensure 
optimal absorption. 

3.2. Space discretization 


QsMj) 


1, if i — j 
0, if i V j 


(19) 


In other words, the <p s j are the usual Lagrange polynomials. The 
total number of nodes in each simplex s £ is equal to D s , so 
that the set {<p s ,i}i=i,D, forms a basis of P s . The total number of 
nodes in the entire domain Q is noted m. 

With the help of these local bases, we can define for each of 
the in nodes i of the domain Q a function <b, such that 


As mentioned earlier, we use the Galerkin finite element 
method [21, 22] to discretize the collective space Q. The main 
reasons for choosing this approach are its capability to manage 
non regular meshes and the possibility to apply h-refinement 
and p-refinement techniques to improve computational effi¬ 
ciency. In this section, we show how to formally derive a linear 
system of equation from the discretization of Eq. (12). 

As customary in quantum mechanics, we note (.[.) the scalar 
product in the space ffttl, C) of complex-valued, square- 
integrable functions, 

(m= f dqf(q)ib(q). (14) 

Jn 


j 4> s ,i(q), if q £ s and i is a node of 5 
q G • VAQ) - | o, otherwise 

( 20 ) 

The functions {0,},=i, m form a basis of our solution space. The 
solution of Eq. (12) can thus be expanded as 

m 

g(q,t) = ^ i g(q i J)(l)i(q). ( 21 ) 

!=1 

Applying the Galerkin finite element method, we search for 
a numerical solution g(q, t) of the form (21) that verifies 

Vf e [0, f max ], Vi £ [1, m\ : Mq)\g(.q, t )) = 0. (22) 
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This process yields a discretized system of m equations with the 
m coefficients giqj, t) as the unknown. It can be written in the 
condensed form 


ihM = [H - ihA]G(t), (23) 

at 

where G(t) denotes the m-dimensional vector of coefficients 
giqj, t) at every node j of the domain Q. The m x m matrices 
M, H and A are defined by 


Mat, = (<t>a(q)\Mq))> 


Aab = {(pa(q)\A(q)<pb(q)), 


H a b = {<t>a(q) I 



+ V(q) 


4>b(q))- 


(24) 


All matrix elements can be computed by applying the basis ex¬ 
pansion (21) to the fields V(q ), B k ,iq) and A(q), that is. 


m 

F(q) = Yj F(q J &(«) with F = V, B kh or A. (25) 

c—l 


The double derivative term in H,j can be integrated by parts 
using the Dirichlet conditions imposed on the boundary d£l. 
The final expression for the matrix elements is 


M a b - 

A a b = 

Hab = 


dq 4> a (q) (/>b(q) 


L 

m ~ 

YMq c ) 

Jn 


dq <f> a (q) <p h iq) </> c (q) 


dq <p a (q) (f> h (q) <pdq) 


rn 

X V(q c ) . 

Jn 

- Yj Yj Bkl(q c) f dq j^-(q)^(q)(f> c (q). 

jri V dn dq k dq, 


(26) 


following prescriptions for the function and its time-derivative, 

8G ^ Gjt + At) - G(t ) a G(t + At) + G(t) 

dt A r ’ 2 ' 1 ; 

Starting from Eq. (23), this numerical scheme yields the fully 
discretized equation 

RxG(t + At) = b{t), (28) 


with 


At At 

R — M + —A + i—H , 
2 2 h 


bit) = 


At At 

M - A - i—H 

2 2 H 


Git). 


(29) 


Using this time discretization scheme, we can show that both 
the norm of the collective wave-function and the average en¬ 
ergy computed from the numerical solution are both constants 
in time if the absorption term is set to zero 


Ilg(?.0ll2 



g*iq,t)giq,t) 


1/2 


cst, 


I 


g\qn) 


H 2 d d 

~2 £ + V( ») 


dq, 


giq,t) = cst. 


(30) 


These properties can be used to test the validity of the numerical 
implementation. 


3.4. Inversion of a linear system 

The collective wave-function at time t is obtained by solving 
the fully-discretized Eq. (28) for the vector Git). This requires 
inverting at each time step a complex, sparse m x m matrices. In 
FELIX these inversions are computed with the iterative QMR 
algorithm without look-ahead as described in [24], The numeri¬ 
cal solution at iteration n is used as the initial guess for iteration 
n + 1. The convergence criterion to stop the iterations for each 
inversion is defined by 


Since the basis functions <p a iq ) are simple polynomials of q. in¬ 
tegrations can be performed analytically, and no quadrature or 
numerical integration scheme is needed. Note that the three ma¬ 
trices obtained here are real, symmetric and sparse. The spar¬ 
sity comes from the fact that the overlap between any two basis 
functions is zero unless at least one element was defined with 
both the corresponding nodes. 

The values of each field at each nodes are inputs of the cal¬ 
culation. To compute the matrix elements, FELIX relies on a 
formal representation of the polynomials. Basis elements <p s j 
are first derived from Eq. (19). This step requires inverting a 
small dense linear system for each simplex. Then derivatives, 
multiplications and integrations of the polynomials involved in 
(Eq. 26) can all be performed formally, so that these operations 
do not generate other errors than those related to the accuracy 
of the polynomial coefficients . 

3.3. Time discretization 

The Crank-Nicolson scheme is used to discretize Eq. (23) in 
time [23]. We recall that the Crank-Nicolson scheme gives the 


\\R x Git + At) - b\\ 2 < e\\b\\ 2 , (31) 

with the tolerance e specified by the user. In order to accelerate 
the inversion, a Jacobi preconditioner is applied to the system 
before the first time iteration. 

Since no look-ahead statement is implemented, the QMR al¬ 
gorithm may occasionally fail to converge at the level of pre¬ 
cision required. In such case, the system is rewritten in the 
(2m x 2m) real form 

/ M + % A -gH \ ( 9Ie(G(f + At)) \ / mbit)) \ 

\ yjH M+^A J\ 3miGit + At)) ) \ Jmibit)) )’ 

(32) 

where 93c and 3m refer to the real and imaginary parts, respec¬ 
tively. This real system is then solved with the Bi-conjugate 
Gradient Stabilized Method as described in [25], 

In practice, the small numerical errors caused by these matrix 
inversions accumulate over time, and can lead to violations of 
the properties (30). However, this numerical error is very small, 
especially if the time span of the time iterations is reasonable 
and the tolerance e is small enough. 
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3.5. Calculation of the flux 

In FELIX, hyper-surfaces are defined as the union of oriented 
faces of an arbitrary list of simplices in the mesh. Note that the 
hyper-surfaces thus defined are not necessarily connected, as 
the simplices need not be adjacent. Given such a hyper-surface 
provided by the user, the code can compute the flux F(f) as 
defined in Eq. (7) through each of the faces f. 

The instantaneous elementary flux / going through an ori¬ 
ented simplex face £' at time t is calculated as 


f (£0 



J k (q,t)dS 


(33) 


where n!- () is the unit vector normal to the simplex face £. J(q. t ) 
is the probability current (6) and N is the dimension of the col¬ 
lective space 1 * * . In order to compute the instantaneous flux, we 
expand the collective wave function g(q, t ) and the inertia tensor 
field Bki(q) on the FE basis using Eqs. (21)-(25). The integral 
in the flux becomes 


f N 

I J k (q. t)dS = ^ B k ,(q u ) 

1=1 u.v.w 

x [9fe(G,(0)3m(G H ,(0) - 3m(G v (f))£Bc(G w (0)] (34) 


P S '. In practice, only two simplices share the interface f. In 
FELIX the value of the integral I^ u ,v,w,i is computed from one 
of these two simplices, which we call the reference simplex. By 
convention, we define the reference simplex as opposite to the 
direction of the normal unit vector, as illustrated in figure 1. 

4. Benchmarks 

In this section, we present a series of benchmark calculations 
that highlight specific features of the code. In each case, the 
analytical solution verifies one or several properties that we use 
to test our numerical implementation. To this purpose, we de¬ 
fine for each case the error between the numerical solution and 
the analytical result, and compute this error as a function of the 
numerical parameters of our calculation, namely 

• the numerical tolerance e for matrix inversions; see 
Eq. (31); 

• the mesh size h, which provides an estimate of the “spa¬ 
tial” resolution of the domain Q; 

• the degree d s of basis polynomials; 

• the time step 5t used in time integration. 


with 

W.w= f Mq)Mq)^r^-ds. (35) 

oqi 

The integral I^ u ,v,w,i is non-zero only if u, v and w are nodes of a 
same simplex containing the edge f. For these non-zero terms, 
the integration is performed formally. This is achieved by first 
computing the polynomial expression of the integrand. Then, 
two successive changes of variables are applied to reduce the 
domain of integration to a simplex of dimension N - 1. Finally, 
we use a trapezoid rule to integrate over time the instantaneous 
flux f(f, t ) to obtain the expression Eq. (7). 



Figure 1: Definition of the reference simplex for a 2-dimensional mesh 


The numerical methods presented here do not enforce the 
continuity of the derivative of the solution at the interface £. 
Therefore, the integral I^ u ,v,w,i ma y not always be defined: for 
two simplices s and s' sharing the interface the value of the 
partial derivative dg/dqi at any point q 6 f may be different 
if computed from the expansion in the basis P s or in the basis 


1 This definition is only valid if the collective space has a dimension strictly 

superior to one. The flux calculation is not enabled in FELIX in the case of 

1-dimensional spaces. 


4.1. Conservation of the Norm 

The conservation of the norm expressed by Eq.(30) is the 
simplest test of our implementation. In this benchmark, the cal¬ 
culation domain G is a 3D cube of size 10 arbitrary units (a.u.). 
The mesh is built by creating a regular grid of equidistant ver¬ 
tices with a mesh size h = 1 a.u.. The position of each vertex is 
then randomly perturbed. A new coordinate q'. is sampled uni¬ 
formly in the interval [q k - ffv, q k + f,h], where is a fluctu¬ 
ation factor set to 15%, and q k is the associated old coordinate. 
Once the vertices are defined, the simplices are determined by 
Delaunay triangulation. In each simplex, a polynomial basis of 
degree d s = 1 is used and the nodes are positioned exactly at 
the vertices. 

Based on this mesh, a free wave packet is propagated during 
a time f max = 15 a.u.. The inertia tensor is diagonal, inversely 
proportional to a mass m = 1.3 a.u. and independent of space. 
The initial wave packet is a 3D isotropic Gaussian centered in 
the middle of the simulation domain and characterized by the 
width cr = fh/ 1.04 a.u.. Finally, we use the time step 5t - 
5.10 4 a.u., and the numerical tolerance for matrix inversions is 
set to its default value of 1CF 15 . 

The error e n0lm of the norm is defined as 


^norm 


lllg(g,0[[2-||g(g,0)|[2l 

lls(ff,0)|| 2 


(36) 


This quantity is computed at different times of the simulation 
and is plotted in figure 2. We note that the error is maintained 
below 10 12 during the whole simulation with this choice for 
the numerical tolerance of matrix inversions. This is consistent 
with the expected property of Eq. (30). As discussed in sec¬ 
tion 3.4 the error <? n0lm comes from the accumulation of errors 
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Figure 2: Error of the norm e norm as a function of time for the free isotropic 
Gaussian wave-packet. 



from matrix inversions at each time step, which explains its in¬ 
crease with time. 

In all subsequent calculations discussed in this paper, we use 
the same numerical tolerance of 10 1 s for matrix inversion and 
the error e norm on the norm will always be below 10 9 . With 
this level of numerical precision, errors coming from matrix 
inversion will always be several orders of magnitude below any 
other source of numerical error examined in this work. 

4.2. Harmonic oscillator potential 

We now turn to the dynamics of a quantum system in an 
isotropic harmonic oscillator (HO) potential in N = 1 and 
N = 2. The advantage of the HO potential is that it provides 
analytical solutions that can be used to test the implementation. 

m at Q fh 

~L3 08 [-20;20] iV 0.15 


Table 1: Characteristics of the ID and 2D harmonic oscillators used in this 
study. 


In the following calculations, the inertia tensor is always di¬ 
agonal, inversely proportional to a mass m and independent of 
space, Bki(q) = 6ki/m. The HO potential being isotropic, it is 
characterized by a single frequency u>, 

V(q) = ^mmq 2 . (37) 


The numerical values adopted for the HO potential are listed 
in table 1. All calculations are performed in a domain Q = 
[-20; 20]^. The same procedure as described in section 4.1 
is used to build an irregular mesh. Dirichlet boundary condi¬ 
tions are enforced at the boundaries <90 of the domain. On the 
other hand, we find that the analytic expressions for the first two 
eigenstates of the HO verify 


VqedQ: 


lg(q)l 

max{ |g(?)|) ?en 


< 2 . 10 -13 . 


(38) 


Therefore, the numerical error coming from the finite size of 
the domain of resolution Q is completely negligible compared 
to the other sources of error under study. 


4.2.1. Ground state of the ID HO 

For any eigenstate of the potential, the modulus | g(q, t)\ of the 
wave function is independent of time. We first test this property 
for the ground state of the ID HO, q = q. The initial wave 
function is taken as 

gm = 0) = exp(-^J, (39) 

where the reduced coordinate Q is defined as 



Calculation is performed up to f max = 32 a.u., which is slightly 
larger than two periods of the complex function g(q,t ). The 
deviation of the modulus of the numerical solution at time t 
from its initial value is measured by the error e mo ±, 

_ ll|G(W)l~|G(f = 0)111- 

emod II |G(f = 0)| |U ’ ' ' 

where the infinity norm |... || LX , for a vector G refers to the max¬ 
imum absolute value of its elements. This error has been com¬ 
puted from a set of calculations with different mesh sizes h and 
time steps 5t. Results are presented in table 2. 


6t\h 

1.0 

0.1 

0.01 

o 

c 

8.089.10~ 4 

1.039.10 " 

7.460.10“ y 

10 -' 

3.780.10 -3 

7.616.10 -6 

1.708.10 6 

1 (T 2 

3.688.10 - 3 

2.603.10- 5 

2.555.10 - 6 

1 (T 3 

3.662.10 -3 

2.626.10 -5 

2.461.10 -6 

1CT 4 

3.661.10 3 

2.619.10 -5 

2.603.10- 6 


Table 2: Error e moc j on the modulus | g(q, r)| for the ground-state of a ID HO. 


We verify that the error e m(ld decreases when refining the cal¬ 
culation in both time and space. For a given mesh size h, the 
error slightly fluctuates with 6t before it finally reaches a con¬ 
verged value. The major part of the error is clearly driven by 
space discretization. 

4.2.2. Sum of two eigenstates for the ID HO 

In this benchmark, the initial state of the system is a sum of 
the first two eigenstates of the ID HO, and can thus be written 

S(e,f = O) = expJ-^J(l + 0. (42) 

The second eigenstate is characterized by a frequency three 
times larger than in the ground-state. Therefore, the motion is 
periodic with the same period as in section 4.2.1. In this case, 
the modulus of the wave function is not constant anymore, since 
the system is not in an eigenstate, but it oscillates between two 
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positions. The full time-dependent solution is given analytically 
by 

g(Q, t ) = exp (-y) (l + Qe~ ,0Jt ). (43) 

We compare the real part of the numerical solution 
93c (g)(q, t ) to its analytic expression using the error defined 
as 

_ l|93e(g num .(f» - 93c(gthe.(0)lloo 
= ii««i, he .(0)ii» • (44) 

This error is computed for several time and space refinements 
and is plotted as a function of time step in figure 3. As in the 
previous benchmark, we observe that the numerical solution 
converges to its analytic value as time and space are refined. 



Figure 3: Error en on the real part of the solution for a periodic motion in a ID 
HO. 


4.2.3. Sum of two eigenstates in 2D 

In this section, the previous study is generalized to the case of 
a 2D HO, q = (q , q'). This benchmark of a 2D case allows us to 
compare the numerical calculation of the flux with its analytic 
expression. The initial wave function reads 

, / Q 2 + Q' 2 \ 

g(Q, 2 . t - 0) - exp I - ——I (1 + Q), (45) 


where Q' is the reduced coordinate associated with q'. Starting 
from this state, the full time-dependent solution reads 


g(Q, Q'J) = exp 


e- + 6- j e - WJt ^ + Ge -K*) (46) 


The system oscillates from one side of the line q — 0 to the 
other side with a period In/co - 31.4 a.u. The derivatives of the 
wave function read 




dg_ 

dq' 


(47) 


mco e 2 +e' 2 .... , ... .a 

e 2 e ■ Q (l + Qe ). 



Figure 4: Definition of the segment [A,B] 


The instantaneous flux of Eq. (33) through an oriented seg¬ 
ment [AB], as depicted in figure 4, is given by 


= V* h sum Sinfojp (Qt sin g_ 0 , t C0S (g)~) 2 
2m 

l yfW L 


f([AB],t) 

x erf (z + Qa cos 6 + Q' A sin 6) J , (48) 

with the error function defined as as usual by 

2 


erf(z) = —— f e ' dx. 
y 7 T Jo 


Frontier 


(49) 



Potential 


1 - 320 . 
1 - 280 . 
1 - 240 . 
— 200 . 
— 160 . 
| 0 
— 80.0 
— 40.0 


Figure 5: Domain of resolution Q for the 2D HO. 

To test our implementation, we first define an arbitrary fron¬ 
tier following an oscillating path around the line q = 0; see fig¬ 
ure 5. The dynamics of the system is computed up to f max = 16 
a.u. At each vertex of the frontier, we calculate the quantities 
93c(g) and <93m (g)/dq at the end of the time evolution. In addi¬ 
tion, our code provides the numerical value of the instantaneous 
flux through each element of the frontier. These three vectors of 
results obtained at the frontier are compared to their respective 
analytic expressions based on the error 


^IronlW) — 


' Hhe.l 


IVthe.llo 


(50) 


These errors have been estimated for a time step 6t = 5.10 4 
and for the following set of different space discretizations, 
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• h = 2,1,0.5,0.2,0.1, with a polynomial basis of degree 
d s = 1. Decreasing the mesh size h is known as li¬ 
re fine merit in the FE approach; 

• h = 2,1,0.5, with a polynomial basis of degree of* = 2. 
For these calculations, each simplex has 6 nodes, of which 
three are positioned at the vertices and the three others in 
the middle of each edge. 

Note that for a same mesh size h, going from linear to quadratic 
polynomials doubles the number of nodes in the mesh. This is 
an example of p-refinement in the FE approach. 



Figure 6: Numerical errors for the real part of the solution, 9fe(g), the spatial 
derivative of the imaginary part of the solution, 53m (g)/dq, and the instanta¬ 
neous flux as a function of the mesh size h for an isotropic 2D HO. 


The results obtained are summarized in figure 6 . Although 
the results are only presented for one value of the time step 5t , 
we obtained similar results with another series of calculations 
using 6t = 1.10 4 . For the three quantities tested, the conver¬ 
gence to their analytic expression is numerically confirmed. We 
note that the convergence is much faster when using polynomi¬ 
als of degree two in the basis. Moreover, the error associated 
with the derivative of the solution is much larger than the one 
associated with the solution itself. In particular, it seems that 
the convergence rate of the flux is limited by the error on the 
derivative. 


4.3. Stokes theorem 

The goal of this benchmark is to test the calculation of the 
flux in a case where no analytic expression of the solution is 
available. This will be achieved by the use of the Stokes theo¬ 
rem of differential geometry. As is well known, if we define an 
enclosed volume V c O, the continuity equation (4) yields the 
following conservation relation 


t’S tokes 


J ^max f* 

dt . 
f=0 JqedV 


J(q, t)dS - 


/1 

JqeV 


\g(q.t)\~dq 


= 0 . 


1=0 


(51) 

We tested this property in the case of a 2-dimensional free wave 
packet. The simulation domain, the mesh, the frontier, and the 
inertia are the same as in section 4.2.3. The initial state is a 


Gaussian function centered at q — (—5; 0) and having a width 
<x = Vfi/1-04 a.u. This function is then multiplied by a plane 
wave of impulsion k = ( 1 ; 0 ) and normalized to one. 

We define the volume V as one half of Q delimited by the 
frontier f and containing the major part of the initial wave 
packet. In this configuration, the Dirichlet condition imposed 
on the boundary of the simulation domain imposes 



t)dS = ^ F(£). 


(52) 


The error estokes on Eq. (51) can therefore be computed at any 
time from the flux and the numerical solution produced by FE¬ 
LIX. 

In this benchmark, we run the simulation up to f max = 3 a.u.. 
At this time, most of the wave packet has crossed the frontier. 
The calculation is performed with a time step 6t = 10 4 and 
repeated for several space discretizations. The results are shown 
in figure 7. The property (51) is verified up to 0.1% of the norm. 



Figure 7: Absolute deviation from the Stokes theorem l^stokesl 


5. Application to the fission of 240 Pu 

To demonstrate the capability of FELIX, we show in this sec¬ 
tion the results obtained by solving the TDGCM+GOA equa¬ 
tion for the neutron-induced fission of a 239 Pu target. In par¬ 
ticular, we emphasize the convergence of such a calculation in 
a realistic setting as well as several features used to increase 
computational efficiency. 

5.1. Description of the calculation 

We solve the TDGCM+GOA equation in the 2D collective 
space spanned by the quadrupole (7/ 20 ) and octupole (q 30 ) mo¬ 
ments of the fissioning system. In the following, qio and qw 
are expressed in fm 2 and fm 3 , respectively. The potential en¬ 
ergy surface Viq) of the compound nucleus 240 Pu in this col¬ 
lective space is obtained by solving the HFB equations with 
a Skyrme energy density in the particle-hole channel and a 
surface-volume, density- dependent pairing energy density; see 
[14] for details. The potential field V(q) is initially computed 


9 






in a domain G 0 characterized by ( 2720 -< 730 ) G [500;60.10 3 ] x 
[-92.10 3 ; +92.10 3 ] (in the units above). In practice, the HFB 
calculation did not converge for every point in the domain; the 
initial grid is thus irregularly spaced and contains 2784 fully 
converged points. At each point, the collective inertia tensor 
B* y' is computed using the perturbative cranking approximation 
of the adiabatic time-dependent Hartree-Fock (ATDHF) theory 
[26], We also determine at each point the expectation value q,\- 
of the Gaussian neck operator, which will be used to determine 
the scission hyper-surface. 

The initial collective state g(q, 0) is defined as the product of 
Gaussian functions in the < 720 - and c/ 3 o-directions. 


g(<720- <730, 0) 



where cr 2 o = 2800 fm 2 and 1 x 30 = 6000 fm 3 are the widths of 
the Gaussian functions, and = 3000 fm 2 , qf^' } = 0 fm 3 
the coordinates of the ground-state. The advantage of such a 
choice is that the initial wave-packet is given analytically and 
does not yield additional numerical errors. 

The resulting wave packet is then multiplied by a plane wave 
characterized by a wave vector k = (2.6.10 - 3 ,0). This last step 
gives the initial state an initial momentum toward positive elon¬ 
gations. It also boosts the average energy of the initial state up 
to roughly 500 keV above the first fission barrier. 

In this work, the isoline q^ = 3.5 mass units defines the 
scission hyper-surface in the collective space. The width of 
the absorption zone is set to w — 8 . 10 3 in the system of units 
adopted. In this area, we impose an average absorption rate of 
r = 20.10 22 s 1 . We solve Eq.(3) up to a time t max = 60.10 22 
s. For the time step values used in table 3 shown in section 
5.3, this corresponds to 60000 and 120000 time iterations. At 
the end of the simulation, the probability for the system to pop¬ 
ulate post-scission configurations is approximately 30%. We 
checked that stopping the time iterations after f max would not 
significantly change the fission fragments yields. 


5.2. Construction of the spatial mesh 

To minimize the computational cost, the time-evolution is not 
performed on a regularly meshed hypercube of the collective 
space. Instead, we use several techniques offered by the finite 
element method to generate a more efficient partition of the do¬ 
main Q. We list below the various steps followed to produce 
the mesh: 

1. Delaunay triangulation - We start with the initial rect¬ 
angular, irregularly spaced domain G 0 mentioned in Sec. 
5.1. We generate a first set of regularly spaced vertices in 
this domain with h = I 120 = 3.1 x /230 the resolution of 
this new mesh Oj. A Delaunay triangulation provides a 
partition of Q[. For every simplex, we choose an interpo¬ 
lation polynomial of degree one. The continuity of each 
field is ensured by placing the three required nodes at the 
vertices. The values of the input fields ( V(q ), Bff, ■ ■ •) are 
then evaluated at each node by linear interpolation in Go. 


2. Absorption areas - The mesh Gj is then extended with 
60 new spatial steps in the lower <720 values, 40 new spatial 
steps in the upper q 20 values and 20 in both <730 directions. 
These extensions define the absorption band. By default, 
the input fields in the absorption band are extrapolated as 
constants based on their values at the edges of the mesh 
G|. The one exception is the potential field in the lower 
<720 region of the absorption band. In this case, we evaluate 
V(q) as a convex parabola continuously connected at <720 = 
500 fm 2 . This simple extrapolation prohibits the system to 
explore oblate shapes during the dynamics. This extended 
mesh is denoted G 2 . 

3. h-refinement and coarsening - We take advantage of the 
flexibility of the finite element method to refine or coarsen 
the mesh G 2 depending on its relevance in the time evo¬ 
lution of the system and the flux calculation. This step 
is needed to improve the numerical precision of the cal¬ 
culation while keeping the computational cost as low as 
possible. After this series of refinement, the new mesh is 
G 3 . 

• In the regions of the domain G 2 where the potential 
takes very large values, the collective wave-function 
g(q, t) will remain very small during the entire time 
evolution. We use the criteria V > Vgs + 35 Mev, 
q N < TO and A H > 170 to automatically detect such 
regions, where we locally coarsen the mesh by dis¬ 
carding two out of three vertices. 

• By contrast, other regions of the domain Qj are im¬ 
portant to the physics of fission, such as near the 
ground-state (where the initial collective state is de¬ 
fined) or near the saddle points. In these regions, we 
apply an additional h-refinement step: each simplex 
is divided in four new elements by adding a vertex in 
the middle of each of its edges. At the boundaries of 
this refinement zone, some simplices may be divided 
differently to ensure the continuity of the fields in the 
new mesh. 

• The most critical region of the potential energy sur¬ 
face is the scission hyper-surface, where one calcu¬ 
lates the total flux. In this area, it is essential to 
maximize the accuracy of the calculation of the func¬ 
tion g(q,t ) and its derivatives. Therefore, we apply 
two successive h- refinement steps near the scission 
hyper-surface. 

• Up to now the domain IT remains rectangular. How¬ 
ever, the areas past and far away the scission hyper¬ 
surface are totally irrelevant to the calculation. We 
thus crop the mesh to retain only the regions of in¬ 
terest. This is done with the following criterion: a 
simplex in the external region is kept only if its dis¬ 
tance to the scission hyper-surface is lower than 10 3 . 

4. p-refinement - We apply one p-refinement step to the 
whole mesh G 3 , that is, after local h-refinement and ca¬ 
reening and cropping. As a result, the degree of interpola¬ 
tion polynomials in each simplex increases to two. This re¬ 
quires three new nodes per simplex, which are positioned 
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in the middle of the simplex edges. The resulting, and final 
mesh, is denoted Q. Starting with h = 397 (in fm units), it 
contains a list L of about 106000 different nodes 
5. Evaluation of fields - Finally, the input fields are again 
evaluated at each node of Q. To do so, results from the 
original calculation on the mesh Clo are linearly interpo¬ 
lated as discussed in the first step. 

The figure 8 shows the final mesh £1 obtained after all the 
aforementioned steps. The main advantage of mesh optimiza¬ 
tion is to reduce significantly the total number of nodes required 
to achieve a given numerical precision. As an example, let us 
consider a two-dimensional, regularly-spaced hypercube mesh 
vV (476 with h = 476 and degree 2 polynomials. We compare 
it with an optimized version Q built from the same initial spa¬ 
tial step. The total number of nodes has been decreased by more 
than 25% even though the mesh is locally up to four times more 
refined in the physically relevant areas. 



Figure 8: Optimized mesh fl obtained with a parameter h = 1190. The red line 
represents the scission hyper-surface, defined arbitrarily by the isoline q n = 
3.5. This figure shows only the region > 0 fm 3 and q 2 o > -10 5 fnT. 


5.3. Results 

We performed a series of calculations with different values 
of the spatial and time resolution parameters h and 6t. The pa¬ 
rameter values h = 793,595,476,397,340,298,264 allow us to 
analyze the convergence in space, while the different time-steps 
6t = 10 3 .5.10 4 (in 10 22 s) control the convergence in time. 

We first consider the total cumulative flux that crosses the 
scission hyper-surface during the whole evolution time. This 
total flux F{t max ) reads 

F(t wax ) — Z F(€, t max ) (54) 


The figure 9 shows the rate of convergence of the flux as a func¬ 
tion of spatial resolution. The relative difference between the 
values obtained for h = 298 and 264 is less than 10 \ We can 
also notice that the calculation is fully converged in time: the 
error on the flux is mostly driven by the spatial resolution. 



Figure 9: Total flux crossing the scission hyper-surface in the interval 
[0; t m ax\ for two different time-steps 5t = 10 3 ,5.10 4 (in 10 -22 s). The 
flux is computed for several meshes characterized by the parameter h = 
793,595.476,397,340,298,264. Note that the x-axis is in log scale. 

Last but not least, we computed the fission fragment mass 
yields as a function of both time and space resolutions. As dis¬ 
cussed in [9], the mass of the fission fragments along the scis¬ 
sion hyper-surface is not necessarily an integer, since both the 
compound nucleus and the fission fragment are described in the 
Hartree-Fock-Bogoliubov approximation, where particle num¬ 
ber is not conserved. As a result, one must take into account 
the fluctuations in particle numbers when comparing the yields 
with experimental data. Following [9], we have convoluted the 
yields coming out of the flux calculation with a Gaussian with a 
width of cr = 3.5 mass units. As customary, the resulting values 
are then normalized to 200. To measure the convergence of the 
yields, we define the quantity <? Y 

e Y = ||7(A) - F ref (A)|U (55) 

The most accurate calculation was obtained with h = 264 and 
6t = 5.10 4 and is chosen as the reference T re f(A). 


6t | h 

793 

595 

476 

397 

340 

298 

264 

io- j 

1.10 

0.69 

0.28 

0.26 

0.13 

0.07 

0.03 

5.10" 4 

1.08 

0.67 

0.27 

0.24 

0.11 

0.05 

0.00 


Table 3: Deviation e\ as a function of time and space resolutions. The param¬ 
eters 6t and h are expressed in 10 22 s and Fermi units respectively. 

We compute the deviation e Y for the different values of our 
numerical parameters. The results are summarized in table 3. 
We first note that the differences on the deviation e Y caused by 
the time resolution are of the order of 0.02. They are much 
smaller than the variations induced by the spatial resolution. 
The evolution of e Y as a function of h shows the convergence 
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in space of our calculation. The values obtained for h = 298 
typically correspond to a numerical error of a few percents for 
masses with Y(A ) >1%. 

6. Program FELIX 

The package FELIX is composed of the following directories 
and files: 

• README: contains detailed instructions to build the solver, 
the tools and their dependencies, and to run the code with 
the examples provided; 

• Makefile: a standard GNU makefile to build the solver 
and the tools; 

• src/: C++ source files of the TDGCM solver and of the 
tools; 

• tools/: additional C++ source files, python and shell 
scripts to handle the inputs and outputs of the TDGCM 
solver; 

• benchmarks/: a few preset inputs and their correspond¬ 
ing outputs; 

• doc/: documentation of the package in DoxyGen format. 

The full Felix package depends on several standard Open 
Source libraries: 

• The TDGCM solver itself requires BLAS, LAPACK, and 
a Fortran compiler with OpenMP support; 

• The documentation requires DoxyGen-1.8.6 or higher; 

• In order to build the full set of tools included in this re¬ 
lease, the user must also install GSL, PETSc, SLEPc and 
Boost. The versions GSL-0.16, PETSc-3.5.2, SLEPc- 
3.5.3 and Boost-1.54 have been used during development. 

6.1. Compilation 

The program is shipped with a Makefile containing a preset 
configuration assuming compilation with the GNU gcc com¬ 
piler on a standard LINUX distribution. The user should expect 
to have to change this Makefile to match his/her own system 
configuration. The Makefile contains some instructions to help 
the user with this configuration step. The different components 
of the package can be compiled separately by typing the fol¬ 
lowing commands: 

• make Doc: generate the DoxyGen documentation in the 
directory, which is located in doc/DoxygenDoc. 

• make Solver: build the TDGCM+GOA solver exe¬ 
cutable. Its name is tdgcmFE and it is located by default 
in src. 

• make Tests: build the executable for the full suite of tests 
included in the package. The name of the executable is 
tdgcmFEtest and is located in tests/src/. 

• make Tools: compile all tools in the directory tools/. 


6.2. Running the solver 

There are two different ways to run the FELIX solver with 
specific input data. 

• If an option file input. opt is available, the user may sim¬ 
ply type 

./src/tdgcmFE input.opt 

• Otherwise, the list of options can directly be passed via the 
command line as: 

./src/tdgcmFE —optionO [valueO] 

—optionl [valuel] . . . 

The available options are discussed below in section 7.2. 

By default, the solver uses every available core on the ma¬ 
chine. As usual, the number of OpenMP threads can be con¬ 
trolled by setting the environment variable OMP_NUM_THREADS. 

6.3. System of units 

By default, the value of the reduced Planck constant h is set 
to 

H = 6.58211928 (l(T 22 MeV.s). (56) 

This imposes a consistency relation between the energy and 
time units that can be used in the code. The most natural choice 
is to set the value of energies in MeV so that time is given in 
units of 10 22 s. This is recommended and is the default setting 
for FELIX. 

If needed, the value of the reduced Planck constant can be 
changed in the file ./src/def ines .h, which allows the user 
to define his/her own set of physical units. Note that the code 
must be entirely re-compiled for such changes to take effect. 
Also, special attention must be paid to setting the units for the 
inertia tensor and coordinates q in a consistent manner. 

7. Inputs and outputs 

7.1. Input files 

FELIX reads most of its input from several ASCII hies. The 
names of these hies begin with the same user-defined prehx, 
and have a specihc extension, as for instance, example. coor, 
example.val, example.geo, etc. The following mandatory 
hies are needed by the solver: 

• example. coor: This hie contains an unsorted list L of 
points in the domain Q defined by their coordinates q (in 
their appropriate unit). In the other inputs hies, a point is 
labeled by its index in the list L. Each line of this hie con¬ 
tains a series of N double numbers separated by a space; 

• example.geo: This hie dehnes the list of all simplices 
of the mesh. Every line contains a series of int integers 
separated by spaces: the hrst integer is the degree of the in¬ 
terpolating polynomial used in the simplex; it is followed 
by the list of the N + 1 vertices of each simplex, followed 
by all the nodes used to dehne the interpolating polynomi¬ 
als. In this version of FELIX the only limitation on the 
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geometry comes from the boundary condition Eq. 11. In¬ 
ternally, this condition is imposed by setting the values of 
the coefficients g(q,, t ) to zero for every node lying on the 
boundary <9fi. The user must therefore design the bound¬ 
ary elements so that the aforementioned property implies 
Eq. 11. 

• example.val: This file contains the values of several 
fields at each point of the list L. The first line starts with 
the special character #, and contains a list of string keys 
that define the fields. Mandatory fields are the potential 
V(q) and the lower part of the inertia tensor. The key for 
the potential is simply v; for the inertia tensor: in 2D, the 
keys are BOO, BIO, Bll; in 3D, they are BOO, BIO, B20, 
B21, B22; etc. Optional fields recognized by the code are 
qN - the expectation value of the Gaussian neck operator, 
and AH - the mass of the heaviest fragment. 

• example. init: It contains the values of the initial wave 
function g(q, t = 0) at each point of the list L. Each line is 
made of two double corresponding to the real and imagi¬ 
nary parts of the wave function. 

The user may also specify additional features and options. 
They are handled through the option list described in section 7.2 
below. Some of these options require one or several additional 
input files: 

• example. front: This file contains a list of oriented 
hyper-surfaces for which the solver will calculate the flux. 
The surfaces are simplex edges and are defined by N ver¬ 
tices. The additional vertex ’’vOpposite” shown in figure 1 
must be specified to set up the orientation. 

• example . matM and example . matH: These two files store 
the sparse matrices M and H defined in Eq. (26). 

• example.opt: A file containing the list of options de¬ 
scribed in section 7.2. 

7.2. Options 

The user can pass a number of options to the solver. Most 
of these options control the numerical parameters of the 
calculation and the frequency at which output data is written 
on disk. 

—help (flag) 

If present, the code only prints a help message and stops the 
execution. Default: Absent. 

—version (flag) 

If present, the code only prints the version number and stops 
the execution. Default: Absent. 

—file (string) 

Prefix name for the input files. All input files must be named as 
prefix.ext where ext is one of the extensions described in 
section 7.1. Default: input. 


—outputDir(string) 

Name of the output directory. If no existing directory is found, 
a new directory is created. Default: . /results 

—dump(integer) 

Number of time iterations between two prints of the solution 
and the flux. Default: 100. 

—refresh(integer) 

Number of time iterations between two displayed lines in the 
standard output. Default: 100. 

—max (integer) 

Maximum number of time iterations for the calculation. 
Default: -1 (no maximum). 

—step(double) 

Time step of the calculation. By default, the unit of the time 
step is 10~ 22 s. Default: 10~ 4 . 

—inversionTol(double) 

Numerical tolerance for matrix inversion; see Eq. (31). Default: 

io- 15 . 

—limit (integer) 

Maximum number of iterations in matrix inversions. Default: 
10 5 . 

—absRate(double) 

Average absorption rate per time unit in the absorption zone. 
For instance, —absRate 20 specifies an absorption rate of 
20.10 22 s -1 , since the basis time unit is 10~ 22 s Default: 0. 

—absWidth(double) 

Width of the absorption zone. The unit depends on the units 
of the collective coordinates q. For any node of the mesh, the 
euclidean distance d to the boundary <90 is computed. A node 
is included in the absorption zone if and only if d <absWidth. 
A negative value will lead to no absorption. Example: Consider 
a rectangular ID space with the axial quadrupole moment given 
in b. Assume the domain is <720 e [0,600] b and absWidth=10. 
Then, all points with <720 £ [0,10] U [590,600] b will be 
included in the absorption band. Default: -1. 

—calcEnergy (flag) 

If this flag is present, the code prints the average energy of the 
solution every dump time iteration. Default: Absent. 

—dumpMat(flag) 

If present, the sparse real matrices M and H are stored in the 
files input. matM and input. matH respectively. Default: 
Absent. 

—readMat(flag) 

If present, the code reads the M and H matrices from the 
files input. matM and input. matH. Note that these matrices 
depend only on the particular mesh of the collective space, but 
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not on time t : for a given mesh, they can be pre-calculated, 
stored on disk using the —dumpMat flag, and re-used in a 
different run. Default: Absent. 

—ISMethod (string=’file’, ’impulsed’ or 
’wavePacket’) 

Method of determination of the initial state g(q, t = 0). 

• If set to ’ f ile ’, the solver reads the initial state from the 
file input. init. 

• If set to ’ impulsed ’, the code reads the initial state con¬ 
tained in the file input. init and multiplies it by a plane 
wave. We note k = kk the wave vector characterizing the 
plane wave, with k the unit vector indicating its direction 
and k the modulus. The user must specify the coordinates 
of k in a file input. k. The modulus k is computed by 
the solver so that the average energy of the initial state 
matches the value provided with the option —ISEnergy. 
The user may provide an initial guess of the parameter k 
via the option —ISLambdaGuess. 

• If set to ’wavePacket’, the solver builds a linear com¬ 
bination of states provided by the user g(q,t — 0) = 
XkVkgk(q)- The user provides the states gkiq) as a set 
of files stated), state_l, etc. A single directory, the 
name of which is set with the option ISStatesDir, must 
contain all the files. The format for the files state_k 
is the same as for input. init. In the expansion of 
the initial state, all weights have a Gaussian dependency 
on the expectation value of the energy of each state: 
ak °c e\p(-(Ek - (E)) 2 /2cr 2 ). The user can tune the pa¬ 
rameter cr via the option ISSigma. The code determines 
the parameter (E) so that the energy of the wave packet 
matches the value provided in the option —ISEnergy. 

Default: ’file’. 

—ISEnergy (double) 

Requested energy of the initial state for the methods 

’impulsed’ and ’wavePacket’. Default: 0. 

—ISSigma (double) 

Gaussian width of the initial state for the method 

’wavePacket’. Default: 1. 

—ISStatesDir (string) 

Directory containing the files state_k required for 

the cconstruction of the initial state with the method 

’wavePacket’. Default: ’ ./eigenStates ’. 

—ISLambdaGuess(double) 

Initial guess for the modulus k of the plane wave multiplying 
the initial state with the method ’impulsed’. If this value is 
negative, the code will initialize k with a default value. Default: 
- 1 . 

—ISNorm (flag) 

If present, the initial state will be normalized to 1 before 


starting time iterations. Default: Absent. 

—frontier(string=’qN’ or ’file’) 

Method of determination of the frontier. If set to ’file’, 
the frontier will be read from the corresponding input file. 
If no such file can be found, the frontier is empty and the 
instantaneous flux is zero. If set to ’qN’, the frontier will be 
computed on the fly as an isoline of the field ’ qN ’. This field 
must then be present as an additional column in the . val file, 
with the key qN. 

Default: ’ f ile ’. 

—frontlso(double) 

Value of the field qN used to define the scission hyper-surface. 
Default: 1. 

—fluxlnst (flag) 

If present, the instantaneous flux through the frontier is 
recorded every dump time iteration. Default: Absent. 

—lumpedMass (flag) 

If present, the lumped mass approximation is applied when 
calculating the M matrix. Default: Absent. 

—pCoeff (double) 

Arbitrary multiplicative factor applied to the potential field 
V(q). Default: 1. 

—bCoeff (double) 

Arbitrary multiplicative factor applied to the inertia tensor field 
Bff. Default: 1. 

7.3. Output files 

All outputs are recorded in a directory that can be specified 
via the option outputDir described in the previous section. 
Upon successful execution of the solver, this directory should 
contain 

• The file example.opt: It contains the list of all options 
used for the run. It could be re-used as an input file for any 
other calculation; 

• The directory gFunction/: Each file g- 

Function.xxx.log in this directory contains the 
solution g after xxx time iterations. These files are for¬ 
matted in the same way as the input file example. init. 
For example, the file gFunction.000000000.log is a 
copy of example. init, unless a renormalization has 
been requested by the user through the use of option 
norm described above. Thanks to this identical format, 
these files can be used in subsequent runs as initial wave- 
functions by simply copying them in place of whatever 
. init file was used. The program FELIX thus has basic 
checkpointing capabilities; 

• The file normDeviation: It contains the deviation 
IlgMIb _ ||g(f = 0)|b as a function of the number of time 
iterations; 
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• The file example.front: This hie contains the frontier 
used for the calculation. It could be either a copy of 
the corresponding input hie or the result produced by the 
solver itself if the frontier is dehned from a qN isoline with 
the option —frontIso=’qN’. 

• The directory flux/: Each hie flux. xxx. log in this di¬ 
rectory contains the time integrated flux F(£) on each ele¬ 
ment £ of the frontier after xxx time iterations. 

Some additional outputs may appear depending on the op¬ 
tions provided: 

• A hie averageEnergy.log: This hie contains the aver¬ 
age energy of the solution as a function of the number of 
time iterations; 

• A hie frontNorm.log: It contains the coordinates of the 
normal vector to each hyper-surface of the frontier; 

• A series of hies fluxlnst .XXX. log in the flux/ direc¬ 
tory: These hies contain the instantaneous flux at the fron¬ 
tier. This flux is not integrated over time. 
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